In [1]:
from numpy import sin, cos, tan, pi, e, exp, log, copy, linspace
from numpy.polynomial.legendre import leggauss
In [2]:
n_list = [2, 4, 8, 10, 20, 30, 50, 100]
In [3]:
f1 = lambda x: exp(x) * cos(x)
f1_sol = -(exp(pi) + 1) / 2
In [4]:
trapezios_simples = lambda f, a, b: (b-a)/2 * (f(a) + f(b))
simpson13_simples = lambda f, a, b: ((b-a)/3) * (f(a) + 4*f((a + b)/2) + f(b))
simpson38_simples = lambda f, a, b: (3/8)*(b-a) * (f(a) + 3*f((2*a + b)/3) + 3*f((a + 2*b)/3) + f(b))
In [5]:
def trapezios_composta(f, a, b, n):
h = (b-a)/n
xi = a
s_int = 0
for i in range(n):
s_int += f(xi) + f(xi+h)
xi += h
s_int *= h/2
return s_int
def simpson13_composta(f, a, b, n):
h = (b-a)/n
x = linspace(a, b, n+1)
s_int = 0
for i in range(0, n, 2):
s_int += f(x[i]) + 4*f(x[i+1]) + f(x[i+2])
s_int *= h/3
return s_int
In [6]:
from sympy import oo # símbolo 'infinito'
def gausslegendre(f, a, b, x_pts, w_pts):
x_gl = copy(x_pts)
w_gl = copy(w_pts)
def gl_sum(f, x_list, w_list):
s_int = 0
for x, w in zip(x_list, w_list):
s_int += w * f(x)
return s_int
if (a == -1 and b == 1): return gl_sum(f, x_gl, w_gl)
elif (a == 0 and b == oo):
x_inf = list(map(lambda x: tan( pi/4 * (1+x)), copy(x_pts)))
w_inf = list(map(lambda w, x: pi/4 * w/(cos(pi/4 * (1+x)))**2, copy(w_pts), copy(x_pts)))
return gl_sum(f, x_inf, w_inf)
else:
h = (b-a)/2
xi = list(map(lambda x: h * (x + 1) + a, x_gl))
return h * gl_sum(f, xi, w_gl)
In [7]:
def erro_rel(est, real):
if real == 0: return abs((est-real)/(est+real)) * 100
else: return abs((est-real)/real) * 100
def aval_simples(f, a, b, real_value):
print('Utilizando os métodos:')
trap_si = trapezios_simples(f, a, b)
print('Trapézio Simples: ' + str(trap_si) + ' Erro Relativo: ' + str(erro_rel(trap_si, real_value)) + ' %')
simps13_si = simpson13_simples(f, a, b)
print('Simpson 1/3 Simples: ' + str(simps13_si) + ' Erro Relativo: ' + str(erro_rel(simps13_si, real_value)) + ' %')
simps38_si = simpson38_simples(f, a, b)
print('Simpson 3/8 Simples: ' + str(simps38_si) + ' Erro Relativo: ' + str(erro_rel(simps38_si, real_value)) + ' %')
In [8]:
def aval_composta(f, a, b, n, x_n, w_n, real_value):
print('Utilizando os métodos: [N = ' + str(n) + '] \n')
trap_c = trapezios_composta(f, a, b, n)
print('Trapézios Composta: ' + str(trap_c) + ' Erro Relativo: ' + str(erro_rel(trap_c, real_value)))
simp_13_c = simpson13_composta(f, a, b, n)
print('Simpson Composta: ' + str(simp_13_c) + ' Erro Relativo: ' + str(erro_rel(simp_13_c, real_value)))
gaule_ab = gausslegendre(f, a, b, x_n, w_n)
print('Gauss-Legendre: ' + str(gaule_ab) + ' Erro Relativo: ' + str(erro_rel(gaule_ab, real_value)))
print('\n')
In [9]:
aval_simples(f1, 0, pi, f1_sol)
In [10]:
for n in n_list:
x_i, w_i = leggauss(n)
aval_composta(f1, 0, pi, n, x_i, w_i, f1_sol)
In [11]:
f2 = lambda x: exp(x)*sin(x) + exp(-x)*cos(x)
f2_sol = lambda x: (exp(x)-exp(-x)) * (sin(x) + cos(x))
x_2 = [0, pi/4, pi/2, 3*pi/4, pi]
h_2 = [0.1, 0.05, 0.01]
In [12]:
df_2pts = lambda f, x, h: (f(x+h) - f(x)) / h
df_3pts = lambda f, x, h: (-f(x + 2*h) + 4*f(x+h) - 3*f(x)) / (2*h)
df_5pts = lambda f, x, h: (-3*f(x+4*h) + 16*f(x+3*h) - 36*f(x+2*h) + 48*f(x+h) - 25*f(x)) / (12*h)
In [13]:
for x in x_2:
print('\nDerivada de f(' + str(x) + ') :' )
d_sol = f2_sol(x)
print('Valor real = ' + str(d_sol))
for h in h_2:
print('com passo \'h\' = ' + str(h) + ' :')
d2r = df_2pts(f2, x, h)
print('Fórmula a 2 pontos: ' + str(d2r) + ' Erro relativo: ' + str(erro_rel(d2r, d_sol)))
d3r = df_3pts(f2, x, h)
print('Fórmula a 3 pontos: ' + str(d3r) + ' Erro relativo: ' + str(erro_rel(d3r, d_sol)))
d5r = df_5pts(f2, x, h)
print('Fórmula a 5 pontos: ' + str(d5r) + ' Erro relativo: ' + str(erro_rel(d5r, d_sol)))
In [14]:
xi, wi = leggauss(100)
In [15]:
f3 = lambda x: x / (1+x)**4
gausslegendre(f3, 0, oo, xi, wi)
Out[15]:
Usando as transformações
$x = \frac{y}{1-y}$
$x = \tan \left[ \frac{\pi}{4} (1+y) \right]$
In [16]:
gausslegendre((lambda y: gausslegendre((lambda x: 1), -(1-y**2)**(1/2), (1-y**2)**(1/2) , xi, wi)), 0, 1, xi, wi)
Out[16]:
In [17]:
gausslegendre((lambda y: gausslegendre(lambda x: e**(-x*y), -(1-y**2)**(1/2), (1-y**2)**(1/2), xi, wi)), 0, 1, xi, wi)
Out[17]: